﻿#include <gtest/gtest.h>

#include "geom/geom_surface_bspline.h"
#include "kernel/spline.h"

#include <algorithm>
#include <iostream>

#include <dbg.h>

using namespace CAGD;

// 测试样条曲面函数
TEST(TestSpline, Surface)
{
    {
        std::vector<std::vector<BPnt>> P(5, std::vector<BPnt>(6));
        for (int i = 0; i < 5; i++)
            for (int j = 0; j < 6; j++)
                P[i][j] = BPnt(i, sin(i + j), j);

        std::vector<double> uknots({0, 0.1, 0.2, 0.4, 0.6, 1, 2, 3});
        std::vector<double> vknots({0, 0.1, 0.15, 0.2, 0.5, 0.8, 1, 2, 3, 3});
        GeomBSplineSurface surface(uknots, vknots, P);

        // 先算两个方向的非零基函数
        int ui = Spline::FindSpan(0.3, 2, uknots);
        std::vector<double> uB = Spline::Basis(ui, 0.3, 2, uknots);

        int vi = Spline::FindSpan(0.6, 3, vknots);
        std::vector<double> vB = Spline::Basis(vi, 0.6, 3, vknots);

        // 组合利用定义计算点
        BPnt S;
        for (int i = 0; i < 3; i++)
        {
            BPnt tmp;
            for (int j = 0; j < 4; j++)
                tmp = tmp + P[ui - 2 + i][vi - 3 + j] * vB[j];
            S = S + tmp * uB[i];
        }

        // 计算曲面上的点，与定义计算的结果比较
        EXPECT_TRUE(surface.Value(0.3, 0.6) == S);
    }

    {
        std::vector<std::vector<BPnt>> P2(5, std::vector<BPnt>(6));
        for (int i = 0; i < 5; i++)
            for (int j = 0; j < 6; j++)
                P2[i][j] = BPnt(i, sin(i + j), j);

        std::vector<double> uknots({0, 0.1, 0.2, 0.4, 0.4, 1, 2, 3});
        std::vector<double> vknots({0, 0.1, 0.1, 0.2, 0.5, 0.8, 1, 2, 3, 3});
        GeomBSplineSurface surface(uknots, vknots, P2);

        double bu, eu, bv, ev;
        surface.Bounds(bu, eu, bv, ev);

        for (int i = 0; i < 100; i++)
        {
            for (int j = 0; j < 100; j++)
            {
                float u = bu + (eu - bu) / 99 * i;
                float v = bv + (ev - bv) / 99 * j;

                int ui = Spline::FindSpan(u, 2, uknots);
                int vi = Spline::FindSpan(v, 3, vknots);

                auto Dp1 = Spline::DerivBasis(2, ui, u, 2, uknots);
                auto Dp2 = Spline::DerivBasis(2, vi, v, 3, vknots);

                // 输入求导后的基函数，计算导数曲面上的点
                auto f = [=](std::vector<double> p1, std::vector<double> p2) -> BPnt {
                    BPnt p;
                    for (int i = 0; i <= 2; i++)
                    {
                        BPnt tmp;
                        for (int j = 0; j <= 3; j++)
                            tmp = tmp + P2[ui - 2 + i][vi - 3 + j] * p2[j];
                        p = p + tmp * p1[i];
                    }
                    return p;
                };

                // 求曲面上的点
                BPnt p = surface.Value(u, v);

                // 求 2 阶导对应的点
                auto P = surface.DerivPoints(2, u, v);

                // 验证所有偏导
                EXPECT_TRUE(P[0][0] == p);
                EXPECT_TRUE(P[0][0] == f(Dp1[0], Dp2[0]));
                EXPECT_TRUE(P[1][0] == f(Dp1[1], Dp2[0]));
                EXPECT_TRUE(P[0][1] == f(Dp1[0], Dp2[1]));
                EXPECT_TRUE(P[1][1] == f(Dp1[1], Dp2[1]));
                EXPECT_TRUE(P[2][0] == f(Dp1[2], Dp2[0]));
                EXPECT_TRUE(P[0][2] == f(Dp1[0], Dp2[2]));
            }
        }
    }
}
